Exact modeling method for discrete finite metamaterial lens 
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We describe an efficient rigorous model suitable for calculating the properties of finite metamate- 
rial samples, which takes into account the discrete structure of metamaterials based on capacitively 
loaded ring resonators. We illustrate how this model applies specifically to a metamaterial lens 
employed in magnetic-resonant imaging. We show that the discrete model reveals the effects which 
can be missed by a continuous model based on effective parameters, and that the results are in close 
agreement with the experimental data. 



I. INTRODUCTION 

For the last decade, metamaterials 0, [|| are in the 
focus of research attention in theoretical and applied 
electrodynamics. Even though no commonly accepted 
definition is available [1, HI, this research direction ex- 
periences a boom encompassing a wide span of areas 
ranging from microwave engineering to nonlinear optics. 
One of the well-known suggestions for applications was 
formulated as a "perfect lens" [jj, making use of neg- 
ative effective material parameters and providing imag- 
ing with subwavelcngth resolution. The idea of super- 
resolution was subsequently analysed and developed in a 
number of ways @, 0, [1], and even realized in practice 
(speaking about three-dimensional systems) using split- 
ring resonators (SRRs) or transmission-line net- 
works [IE [3. 

Arguably the closest approach to practice offered by 
metamaterials, is related to magnetic resonance imaging 
(MRI). For example, rotational resonance of magnetoin- 
ductive waves jl3| was suggested for parametric ampli- 
fication of MRI signals 14J or enhanced detection with 
flexible ring resonators 15| . Alternatively, applications 
based on 'swiss- rolls' [161 or wire media [17| channelling 



were put forward [16J, [18( . Naturally, superlens concept 
is also promising in this area: specifically for MRI, an 
isotropic lens based on capacitively loaded single ring res- 
onators was designed and experimentally tested [10| . 

Such a metamaterial lens is intended to operate at 
the value of effective permeability fj, = —1. In theory, 
for modelling such metamaterials (based on SRRs), one 
can exploit an effective medium approach, taking care 
of numerous limitations related to general restrictions of 
homogenization [lj| as well as to specific peculiarities 
caused by resonant nature of the structural elements [2(| • 
Universally, all the structure details (size of the elements 
and lattice constants) must be much smaller than the 
wavelength; while the total number of elements in meta- 
material should be sufficiently large to make homogeniza- 
tion meaningful. In addition, spatial dispersion effects 
can be rather remarkable in metamaterials, and impose 
further restrictions on effective medium treatment, pro- 
hibiting that in certain frequency bands (20j . 

Unless one opts for a completely numerical homoge- 



nization method 21[, generally applicable to almost ar- 
bitrary structures, a model have to be developed to de- 
scribe adequately the effective medium properties. Quite 
general approach (20| for homogenization of resonant 
metamaterials can be applied to a variety of metama- 
terials including those which combine different clement 
types. However, this relies on the dipole approxima- 
tion for the interaction of elements in the lattice, which 
may not be always valid. For example, mutual inter- 
action of the circular currents close to each other sig- 
nificantly differs from dipole interaction, which becomes 
relevant for dense metamaterials. The first rigorous anal- 
ysis of such metamaterials was given early in Ref. [22| . 
where the effective permeability has been derived given 
the properties of individual elements and lattice param- 
eters, through the classical procedure of averaging the 
microscopic fields. In that approach, mutual inductances 
between a large number of neighbours are taken into ac- 
count, revealing the importance of lattice effects. This 
approach, however, is limited to quasi-static conditions, 
and requires wavelength to be much larger than any 
structural details. Later on, a rigorous method was elab- 
orated for isotropic lattices of resonant rings [23[ , which 
accounts for spatial dispersion. On the other hand, the 
latter approach employs a nearest neighbour approxima- 
tion as otherwise full analytical treatment becomes im- 
practical. 

The above theoretical methods provide the effective 
parameters for "infinite" structures (which in practice 
implies the structures sufficiently large in all three di- 
mensions). The lens of Ref. [ldj . however, contains just 
a few elements across the slab. Specifically for this case, 
a method was developed to calculate the transmission 
properties for a thin infinite slab (24j; furthermore, it was 
shown that similar results can be obtain by considering 
an equivalent slab with the effective medium parameters 
as obtained in Ref. [23| . 

Nevertheless, it is clear that a number of peculiar ef- 
fects caused by the discrete structure of the lens as well 
as its finite size, cannot be reliably assessed with the 
above models, as the lens is too small for an effective 
medium treatment. On the other hand, it is large enough 
to make an analysis with full-wave commercial software 
practically impossible. For this reason, here we develop a 
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finite model to calculate lens properties, which explicitly 
takes all the structural details into account. The goal of 
this paper is to describe this modelling approach in de- 
tail, and to illustrate that indeed it does reveal some fea- 
tures which are missed by the continuous modelling. We 
should note, though, that while the model is described 
in connection to one particular structure, the approach 
applied here is generally applicable to any realistic SRR- 
based metamaterial, and therefore is useful for a wide 
range of applications. 



II. GEOMETRY OF THE PROBLEM 

The metamaterial lens described in Ref. [l(| is com- 
posed of capacitively loaded rings (CLRs) periodically 
arranged in an isotropic three-dimensional lattice with 
the lattice constant a = 1.5 cm. The lens features three 
planes of 18 by 18 CLRs interlayered with orthogonal 
segments providing two mutually orthogonal sets of two 
layers 17 by 18 CLRs each (see Fig.rjjfor clarity), which 
makes it up to roughly 2200 CLRs. This lens can be op- 
tionally extended by an extra 3D-laycr, resulting in hav- 
ing four 18 x 18 layers interlaced with the two orthogonal 







FIG. 1: Photograph of the quasi-magnetostatic metamaterial 
lens analysed in this paper. 



FIG. 2: (a) Sketch of the CLR resonator; (b) Scheme of the 
lens with the corresponding coordinate system. 



subsystems of 3 by 17 by 18 CLRs, amounting to about 
3130 elements. Overall dimensions of the (non-extended) 
lens are thus 27 x 27 x 3 cm. 

The CLRs themselves (Fig. [2^,) are made of copper 
through etching metallic strips on a dielectric board. The 
mean radius ro of the CLRs is 0.49 cm (r^/a = 0.33) and 
the strip width w is 0.22cm (w/a = 0.15). The CLRs 
are loaded with lumped non-magnetic 470 pF capacitors. 
The sclf-inductancc of the CLRs, L = Wq/C = 13.5 nH, 
has been obtained from the measured value of the fre- 
quency of resonance in free space, equal to 63.28 MHz 
(koa = 0.02). By measurement of the quality factor 
of the resonator the resistance has been estimated as 
R = 0.0465 Ohm, which includes the effects of both the 
ring and the capacitor. 

We reserve the standard coordinate system (x, y, z) 
for discussions on the level of geometry of one ring and 
their mutual interactions, as relevant for the next sec- 
tion. When referring the overall lens geometry, we define 
supplementary coordinate system (X, Y, Z) so that the 
lens geometrical centre is placed at the coordinate origin, 
and the Y axis is perpendicular to the lens as slab ( "lens 
axis"), while the long edges of the lens are parallel to X 
and Z axes (Fig. [2Jd) . Note that the coordinate origin is 
between the rings, and thus the lens is completely sym- 
metric with respect to the coordinate origin, all the axes 
and all the coordinate planes (in the analysis, we neglect 
minor asymmetry occurring in the real lens caused by 
specific assembly details, e.g. resulting from substrate 
thickness, as these deviations are of the same order as 
unavoidable production inaccuracy). For the three mu- 
tually orthogonal sets of CLRs, we will refer to as X-rings, 
Y-rings or Z-rings, depending on whether the rings' nor- 
mals are along X, Y or Z axes. The lens, therefore, con- 
tains 612 of either X-rings and Z-rings, and 972 Y-rings. 
We will also introduce a consecutive numbering of all the 
CLRs with a single index. The so called input and output 
surfaces of the lens correspond to Y = =pl.5 cm, while the 
theoretical source and image planes are at Y — =p3.0cm. 



3 



III. THEORETICAL MODEL 

For the analysis of the lens response to the external 
field, we consider an ideal cubic lattice of L-C circuits 
supporting current. With the time convention as I oc 
exp(ju;i), each of the currents is governed by equation 



ZqIu 



(1) 



where the self-impedance Zo = (R + juiL + 1/ (jwC)) is 
determined by the resistance R, self-inductance L and 
sclf-capacitancc C of the single CLR, while <&„ represents 
the total magnetic flux through the considered ring which 
can be written as 



<r> — * oxt J- \^ <j> — <j> cxt 4. \^ m T 



(2) 



where $° xt is the magnetic flux from external sources 
and M nm are the mutual inductances between the rings 
n and m. Combining Eq. fl} and Eq. ((2]) we obtain 



Z I 



(3) 



with Z„ 



Zn . Z„ 



ju)M nm , which is a system of 



linear equations for unknown currents, provided that the 
external sources are known. 

Mutual inductance between the flat rings (which is the 
case under consideration) carrying the currents /„ and I m 
uniform along the ring contour is, most generally [25| . 



M„ 



K„(r)-K m (r') 



dSdS' (4) 
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where K represent surface current densities; we assume 
that these follow Maxwcllian distribution across the strip, 
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(5) 

Clearly, such integration is not ideally suited for nu- 
merical calculation. In the first approximation, mutual 
inductance between CLRs can be estimated with the one 
between linear currents (double linear integration along 
the equivalent ring contour), but for close CLRs this does 
not give a good precision. 

However, a trick is that the result of surface integration 
according to Eq. @ can be approximated with a good 
precision through an average mutual inductance between 
two pairs of circular currents [26| . This way, each flat ring 
can be represented by a pair of coaxial circular currents 
of radii ro ± 7^/2, and the sought mutual inductance is 
calculated as an average between the four corresponding 
linear ones: 



Mnm — (i^m + ^nm + ^nm + ^nm) 1^1 (6) 




FIG. 3: Geometry of the linear currents for parallel or or- 
thogonal ring orientations, as relevant for mutual inductance 
calculations. 



which essentially decreases calculation time. The value 
of particular parameter 7 depends on the ring geometry, 
but does not depend remarkably on the relative orienta- 
tion and distance between the CLRs (within the limits 
of lens structure) . For the particular parameters consid- 
ered here, 7 « 0.7 was numerically found to give a good 
match to the precise integration ^ (while 7 = 1 would 
correspond to the edges of the strip). 

To achieve a faster calculation for the "linear" mutual 
inductance itself, we further note that it can be easily 
evaluated [25| by integrating the vector potential A along 
the current contour I: 

A„-dl m , (7) 



L, 



In 



with vector potential itself having only A v component 
(assuming that the source ring is placed at the coordinate 
origin with the normal along z axis, sec Fig. [3]), which 
can be obtained with the help of elliptic integrals as [25[ 




2 - k 2 ) K{x) - 2£{x) 



(8) 



with 



Ar a p 



(p + r ) 2 + 2 s 



being the argument of complete elliptic integrals 



7T/2 



it/2 



£= Vi" 



x 2 sin 2 6 d0, 



x 2 sin 2 9 



Given the fact that the fast pre-defined routines for ellip- 
tic integrals are available in a number of computational 
platforms (e.g. Matlab®), this effectively reduces double 
integration to a single one. Thus, finally, only four linear 
integrals like ([7]) arc required to approximate the exact 
value of a double integration like (H}. 
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IV. NUMERICAL IMPLEMENTATION 

To analyse the response of the lens to an external field 
source, the key step lies in solving system ((3|). To do so, 
we need to know the matrix of mutual inductances M. 
This matrix is only determined by the geometry of the 
rings arrangement inside the lens and can be calculated 
once for a given lens geometry, while the impedance ma- 
trix Z can be then obtained for all frequencies as shown 
after ©. 

For such a lens as described above, having 2196 rings, 
the matrix contains almost 5 million values, and filling 
those with a direct calculation would be rather time- 
consuming even with a simplified integration described 
in the previous section. However, obvious reciprocity 
(M nm = M mn ) and symmetry properties of the lens al- 
low for a great simplification of matrix filling. Indeed, 
the lens is symmetric with respect to J, 7 and Z axes 
as well as to XY, YZ and ZX planes. This implies, in 
particular, that the mutual inductances between X-rings 
and Y-rings are all the same as between Z-rings and Y- 
rings. Furthermore, as all the rings are identical, induc- 
tance between them is only determined by their mutual 
orientation and spatial offsets Ax, Ay, Az (see Fig. [3j) , 
and for parallel rings even Ax is equivalent to Ay. Ex- 
plicitly, integration for the mutual inductances between 
the parallel rings is performed according to 



L P (A6,A,) = / A / 2(r2 + A6c ° SQ) d 



Ap 



(9) 



Ap 



>c 2 = 



{Ab) 2 + 2r 2 A&cosa, 
4riAp 



where A6 = ^(Az) 2 + (Ay) 2 
orthogonal mutual orientation 



(Ap + n ) 2 + (Az) 2 

and for the rings with 



2tt 



L°(Ax,Ay,Az) 



r 2 Ay cos a 

A v da, 

Ap 



(10) 



Ap = sj (Ax — r 2 sina) 2 + (Ay) 2 , 
2 4nAp 

K — . 

(Ap + ri) 2 + (Az — r 2 cos a) 2 

In the above equations, we imply a general case that the 
radii of the two rings (r\ and r 2 ) can be different. 

Thus, a number of ring pairs within the lens share 
the same value of mutual inductance, so it is only neces- 
sary to calculate a full set of non-equivalent mutuals and 
then assign those values depending on the mutual offsets. 
With the particular lens considered here, there are only 
1924 independent inductances for the parallel ring ori- 
entation, and 1668 for the orthogonal one, so the total 
number of calculations © is 3592 — orders of magnitude 
smaller than the number of matrix elements. This way, 



the entire matrix can be filled in a matter of seconds on 
an ordinary PC. 

Another preliminary step is to determine the external 
flux $ cxt imposed to each ring by a given source. For a 
homogeneous field or a plane wave excitation, calculation 
is straightforward with the known coordinates of each 
ring: 



<PT = nr 2 n B„ ■ n, 



(11) 



where n is a ring normal while magnetic field B„ can be 
evaluated at the ring centre as the field variation across 
the ring is negligible. 

In practice, the lens is typically used along with exci- 
tation / measuring coils employed in MRI practice. In 
that case, instead of calculating the field produced by a 
coil over each ring (which is further complicated as this 
field is not uniform across the ring), it is much easier to 
obtain the flux directly 
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FIG. 4: Frequency dependence of the real (a) and imaginary 
(b) parts of the impedance measured by a 3-inch coil placed 
at the image plane (Y = —3 cm). 
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in terms of mutual inductance M% between the coil and 
each ring, which can be calculated with the same method 
as the one between the rings. Above, I c = Zn+x is the 
total current induced in the coil by the external voltage 
source as well as by the lens. Imposing a given voltage 
V c to the coil with the self- impedance Z c , we can include 
the coil mutual impedances into system ((3j> , modified as 



Z 1 = V 

\}L0M m 



(13) 



for 
for 



1 < n < N 

71 = N +1 



with V n = V c S„, t N and N being the total number of rings 
in the lens. Clearly, additional coils, if necessary, can be 
included by extending the matrix system in an analogous 
way. 

After the above procedures, it is finally possible to 
solve the systems ([3]) or (|T3| obtaining currents /„ in 
each ring for any given excitation. With these known, 
it is further possible to calculate any desired response of 
the lens, such as magnetic field produced by the lens (us- 
ing standard Biot-Savart expressions) or impedance as 
measured by the MRI coil, 



7COil 



In 
Z c 



(14) 



V. RESULTS AND DISCUSSION 

Armed with the above precise method, we can have 
a detailed look into lens features and response to vari- 
ous external field sources. In previous work [24| it was 
concluded that the accurate model, developed for an 2D- 
infinitc slab with the same structure and thickness as the 
real lens, is capable of predicting the observations made 
in connection to lens use in MRI practice. In a typical 
setup, a coil of 3 inch in diameter is placed parallel to the 
lens interface at the source plane, Y = —3 cm (that is, 
at a distance 1.5 cm, equal to one half of the lens thick- 
ness from the lens surface). The super- lens behaviour 
implies that the magnetic field produced by the coil, is 
then reproduced in the space behind the image plane 
(Y = 3 cm) , as if the coil itself were present in place of 
the image. 

A straightforward example to test the developed model 
and to compare it with practice as well as with earlier 
models, is to evaluate the impedance as measured by a 
coil in front of the lens, depending on frequency. In the 
discrete model, it is given by Eq. (TH| . while with the 
continuous model (for an infinite homogeneous slab with 
an appropriate effective permeability [23]) it can be nu- 
merically calculated as 



Z coil = -— Re / E 1 dl 



< 




(15) 



distance (cm) 

FIG. 5: Axial component Hy of the total magnetic field ob- 
served behind the lens surface along the lens axis (1) or along 
the parallel line (2) slightly displaced in X and Z direction so 
that it passes through the centre of one ring (see the inset for 
the labels of the axes). Comparison between the two models 
when the lens is excited by 3-inch coil, centred with respect 
to the lens axis and placed at Y = —1.5 cm. 



where E r is electric field reflected by the lens. The two 
modelling results are compared in Fig. 0] along with the 
experimental data. Although there is no exact quantita- 
tive matching to the measured data, it is clear that the 
frequency dependence provided by the discrete model is 
closer to experiment than that of the continuous calcula- 
tion. On the other hand, we can conclude that the latter 
already provides qualitatively suitable picture, predict- 
ing an overall pattern of the impedance frequency de- 
pendence. 

With both the continuous model [24| and the model 
developed above, it is easy to obtain the axial magnetic 
field Hy behind the lens for a given excitation. Compar- 
ison between the predictions of the two models is shown 
in Fig. One can see that at distances smaller than 
about one lattice constant (a = 1.5 cm), Hy is essen- 
tially inhomogeneous as the near-field of the individual 
rings dominates, so that the total field is quite different 
whether traced along the lens axis (which passes between 
the rings) or along a line that passes through a ring cen- 
tre, while both are remarkably different from the continu- 
ous model. This is an obvious consequence of the discrete 
lens structure, which cannot be revealed by a homoge- 
nized model but is apparent in practice. At distances 
larger than approximately one lattice constant (1.5cm), 
the field observed along the two axes converge, and are 
qualitatively similar to the continuous model with a fair 
numerical agreement (see Fig. [5]). 

Another peculiarity arising from the discrete structure 
is related to the spatial resolution of the lens in the X-Z 
plane. Evidently, a lens cannot resolve any details which 
are separated by distances of the order of lattice constant. 
To identify the actual limitation, we test the magnetic 
field distributions originating from using the coils of var- 
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FIG. 6: Axial component Hy of the magnetic field observed 
behind the lens surface. Horizontal axis corresponds to the 
lens surface (parallel to X-Z plane), while the vertical one is 
normal to the lens (Y). Only one half of the symmetrical field 
spatial distribution is presented; normalized magnitudes are 
shown in logarithmic scale between 10 (white) and 0.1 (black) 
A/m. Excitation with coils of different radii (0.5 cm, 1.5 cm, 
3 cm and 4.5 cm), centred with respect to the lens axis and 
positioned at Y = —1.5 cm. 



assess the source location and size. A reasonable picture 
is obtained for a 4.5 cm coil radius, where the field farther 
than the image plane looks as expected with super-lens 
performance (Fig. We can therefore conclude that 
spatial resolution of the discrete lens can be assumed to 
be of the order of 5 lattice constants. This observation is 
in good agreement with the general concerns regarding 
the lattice effects in metamaterials (22j . 

With the above examples, we clearly demonstrate that 
the exact model described in this paper, is suitable for a 
reliable description of the metamaterial lens, and makes 
it possible to predict specific observations which might 
be missed by a continuous model. 

Certainly, the above methodology is not restricted to 
the particular lens geometry and can be perfectly used for 
any metamaterials designed with CLRs or SRRs, whether 
isotropic or anisotropic, and also arbitrarily small in size. 
The only limitation is that for very large number of el- 
ements, numerical evaluation on conventional computers 
may fail, specifically as far as allocating space for huge 
impedance matrices, and inverting these, is concerned. 
However, in metamaterial research it rarely comes to 
samples that large, and, on the other hand, when it 
comes, then there are good reasons to expect that con- 
tinuous models will work sufficiently fine. 

In contrary, for small metamaterials typically consid- 
ered for practical use, modelling this way provides an 
invaluable insight into their properties and leads to reli- 
able predictions. 



ious small radii (Fig. [6]). For excitations with a coil of the 
ring size, the entire lens is dominated by standing mag- 
netoinductive waves [27j], and the field pattern does not 
suggest any hints for resolving the source (Fig. [TjVi). In- 
deed, practically the same field pattern is observed with a 
three times larger coil, two lattice constants in diameter 
(Fig. [BJd). With a still larger coil, encompassing four lat- 
tice constants, one may argue that the pattern starts to 
clarify (Fig. Et), though still it cannot be reliably used to 
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